* ----------- *
* Making Data *
* ----------- *

clear all

set scheme lean2

set seed 123

set obs 4500

gen x = _n

gen year = ceil(x / 45)

gen period = year > 60

gen time = year - 60

gen r1 = runiform()

cap drop singular
gen singular = 0

forvalues i = 1 / 60 {
	replace singular = 1 if r1 < `i' / 200 & year == `i'
}


forvalues i = 61 / 100 {
	replace singular = 1 if r1 < .3 & year == `i'
}


gen r2 = runiform() if singular == 0


forvalues t = 0 / 40 {
	replace singular = 1 if r2 < `t' / 40 & time == `t'
}


* -------------------------------------- *
* Regression Results With Yearly Cutoffs *
* -------------------------------------- *

matrix define A = J(61, 3, 0)
*matrix define B = J(61, 3, 0)
*matrix define C = J(61, 3, 0)

forvalues t = 20 / 80 {
	
	cap drop time
	cap drop period
	
	gen time = year - `t'
	gen period = t > 0
	
	qui reg singular c.time##i.period, cluster(time)
	
	local row = `t' - 19
	
	*matrix C[`row',1] = `t'
	*matrix C[`row',2] = `e(r2_a)'
	*matrix C[`row',3] = `e(rmse)'

	
	qui lincom 1.period#c.time
	matrix A[`row',1] = `t'
	matrix A[`row',2] = `r(estimate)'
	matrix A[`row',3] = `r(p)'
	
	*qui lincom 1.period
	*matrix B[`row',1] = `t'
	*matrix B[`row',2] = `r(estimate)'
	*matrix B[`row',3] = `r(p)'
}

cap drop A1 A2 A3
svmat A
twoway ///
	(scatter A2 A1 if A3 < 0.05, msym(o)) ///
	(scatter A2 A1 if A3 > 0.05, msym(oh)) ///
	, legend(off) ///
	xlab(20 "1820" 40 "1840" 60 "1860" 80 "1880") ///
	xtitle("Cutoff Year") ///
	ytitle("Slope Change")
graph export "${output}FigS3_simulation_slope.pdf", as(pdf) replace





* ------------------------- *
* Graphs of What's Going On *
* ------------------------- *

collapse singular, by(year)

twoway ///
	(lfitci singular year if year < 61, lpattern(solid)) ///
	(lfitci singular year if year >= 61, lpattern(solid)) ///
	(scatter singular year, msym(o)) ///
	, legend(off) ///
	xlab(0 "1800" 20 "1820" 40 "1840" 60 "1860" 80 "1880" 100 "1900", labsize(vlarge)) ///
	ylab(, labsize(vlarge)) ///
	xline(1860) ///
	xtitle("Year", size(vlarge)) ///
	title("1860 Cutoff", size(huge)) ///
	name(c1860, replace) nodraw
*graph export "${output}simulation_1860.png", as(png) replace
	

twoway ///
	(lfitci singular year if year < 31, lpattern(solid)) ///
	(lfitci singular year if year >= 31, lpattern(solid)) ///
	(scatter singular year, msym(o)) ///
	, legend(off) ///
	xlab(0 "1800" 20 "1820" 40 "1840" 60 "1860" 80 "1880" 100 "1900", labsize(vlarge)) ///
	ylab(, labsize(vlarge)) ///
	xline(1830) ///
	xtitle("Year", size(vlarge)) ///
	title("1830 Cutoff", size(huge)) ///
	name(c1830, replace) nodraw
*graph export "${output}simulation_1830.png", as(png) replace
		
	
twoway ///
	(lfitci singular year if year < 80, lpattern(solid)) ///
	(lfitci singular year if year >= 80, lpattern(solid)) ///
	(scatter singular year, msym(o)) ///
	, legend(off) ///
	xlab(0 "1800" 20 "1820" 40 "1840" 60 "1860" 80 "1880" 100 "1900", labsize(vlarge)) ///
	ylab(, labsize(vlarge)) ///
	xtitle("Year", size(vlarge)) ///
	title("1880 Cutoff", size(huge)) ///
	name(c1880, replace) nodraw
*graph export "${output}simulation_1880.png", as(png) replace		
	
graph combine c1860 c1830 c1880, row(1) xsize(3) ysize(1)
graph export "${output}FigS2_simulation_cutpoints.pdf", as(pdf) replace		




